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ABSTRACT 

The motion of a black hole about the centre of gravity of its host galaxy induces 
a strong response from the surrounding stellar population. We treat the case of a 
harmonic potential analytically and show that half of the stars on circular orbits in 
that potential shift to an orbit of lower energy, while the other half receive a positive 
boost and recede to a larger radius. The black hole itself remains on an orbit of fixed 
amplitude and merely acts as a catalyst of evolution of the stellar energy distribution 
function f{E). We show that this effect is operative out to a radius of w 3 to 4 times 
the hole's influence radius, Rhh- We use numerical integration to explore more fully 
the response of a stellar distribution to black hole motion. We consider orbits in a 
logarithmic potential and compare the response of stars on circular orbits, to the 
situation of a 'warm' and 'hot' (isotropic) stellar velocity field. While features seen in 
density maps are now wiped out, the kinematic signature of black hole motion still 
imprints the stellar line-of-sight mean velocity to a magnitude ~ 18% the local root 
mean-square velocity dispersion a. 

Key words: numerical method: N-body; galaxies, gravitational dynamics 



1 Introduction 

Black hole dynamics in galactic nuclei has attracted much 
attention for many years (e.g., Begelman et al. 1984; Kor- 
mendy & Richstone 1995; Merritt 2006 for a recent review). 
The influence of a black hole on its surrounding stars is felt 
first through the large velocity dispersion and rapid orbital 
motion of the inner-most cluster stars {a ~ fid ;$ 10^ km/s). 
This sets a scale < GMbh/cr^ 0.015 - 0.019 pc for the 
Milky Way, henceforth MW) within which high-angle scat- 
tering or stellar stripping and disruption may take place. 
For the MW, low-impact parameter star-BH encounters are 
likely given the high density of p ~ 10^ Mq /pc^ within a 
radius of ~ 10 pc (see e.g. Yu & Tremaine 2003; O'Leary 
& Loeb 2006; see also Freitag et al. 2006 for a numerical 
approach to this phenomenon). Star-BH scattering occur- 
ring over a relaxation time (Preto et al. 2005 and references 
therein; Binney & Tremaine 1987) leads to the formation 
of a Bahcall-Wolf stellar cusp of density p* ~ where 7 
falls in the range 3/2 to 7/4 (Bahcall & Wolf 1977). Genzel 
et al. (2003) modeled the kinematics of the inner few par- 
sees of Sgr A* with a mass profile p* ~ t""^ *, suggestive of 
a strong interplay between the black hole and the central 
stellar cusp. More recently, Schodel et al. (2007) presented 
a double power-law fit to the data, where the power index 
~ 1.2 inside a breaking radius rtr, and ~ 7/4 outside, where 



ri,r — 0.2 pc. This is indicative of on-going evolution inside 
Vbr not accounted for in the Bahcall-Wolf solution. 

Most, if not all, studies of galactic nuclei dynamics as- 
sume a fixed black hole (or black hole binary) at the centre 
of coordinates. Genzel et al. (1997) had set a constraint 
of < 10 km/s for the speed of the black hole relatively to 
the galactic plane, a constraint later refined to < 2 km/s 
(Backer & Sramek 1999; Reid & Brunthaler 2004). SteUar 
dynamics on scales of ~ few pc surrounding Sgr A* is 
complex however, and the angular momentum distribution 
on that scale is a prime example of this complexity (Genzel 
et al. 2003). Reid et al. (2006) used maser emission maps to 
compute the mean velocity of 15 SiO emitters relatively to 
Sgr A*. They compute a mean (three-dimensional) velocity 
of up to 45 km/s, a result obtained from sampling a volume 
of ~ 1 pc about the centrcj^ This raises the possibility that 
stars within the central stellar cusp experience significant 
streaming motion with respect to Sgr A*. The breaking 
radius rt,- ^ 0.2 pc is suggestive of uncertain dynamics 
on that scale. Random, 'Brownian' black hole motion may 
result from the expected high-defiection angle encounters 
(Merritt 2005; Merritt et al. 2007). Here we take another 

^ Statistical root-n noise ~ 25% remains large owing to the small 
number of sources but is inconsequent to the argument being 
developed here. 
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approach, and ask what net effect a black hole set on a 
regular orbit will have on the stars. In doing so, we aim to 
fill an apparent gap in the modeling of black hole dynamics 
in dense nuclei, by relaxing further the constraint that the 
hole be held fixed at the centre of coordinates. 



A rough calculation will help to get some orientation 
into the problem. Consider a point mass falling from rest 
from a radius Ro in the background potential of the MW 
stellar cusp. Let the radial mass profile of the cusp p*(r) oc 
r~^^^, consistent with MW kinematic data. If we define the 
black hole radius of infiuence ~ 1 pc to be the radius where 
the integrated mass M*(< r) = the black hole mass ~ 3 to 
4 X 10® Mq (Genzel et al. 2003; Ghez et al. 2005), then Ro 
may be expressed in terms of the maximum black hole speed 
in the MW potential as 



max{ii} 
lOOkm/s 



Ro 
1 pc 



For a maximum velocity in the range 10 to 40 km/s we find 
Ro — 0.3 — 0.5 pc, or the same fraction of its radius of in- 
fiuenc^ We ask what impact this motion might have on 
the surrounding stars. To proceed further, let us focus on a 
circular stellar orbit outside Ro in the combined potential 
of the black hole and an axisymmetric galaxy. Fig. [l] gives 
a clue to the analysis. When the black hole is at rest at the 
centre of coordinates, the star continues on a closed circu- 
lar orbit of radius r and constant velocity v. We now set 
the black hole on a radial path of amplitude Ro down the 
horizontal x-axis. Without loss of generality, let the angular 
frequency of the stellar orbit be to' , and that of the black 
hole ui ^ uj' . The ratio uj/lu' 1 is otherwise unbounded. 
The net force F acting on the star can always be expressed 
as the sum of a radial component Fr and a force parallel to 
the X-axis which we take to be of the form Fx cos{ujt + ip); 
clearly the constant Fx — when Ro — 0. The net me- 
chanical work done on the star by the black hole as the star 
completes one orbit is computed from the integral 

5W = j F-vdt = J FxVsm{u't)cos{Lot + •p)dt (1) 

where if is the relative phase between the stellar and black 
hole orbits. The result of integrating ([T]) is best set in terms 
of the variable v = lj/uj' as 



5W 

rFx 



- 1 



[cos(27rz^ + tfi) — cos{lp)] 



when u > 1, and 



SW 

rFx 



— 77 sin((^) 



(2) 



(3) 



when = 1. It is a simple exercise to show that this last 
expression is recovered from ^ in the limit !^ — > 1"*". Equa- 
tion (3) embodies the essential feature, which is that SW 



These figures are robust to details of the stellar cusp mass pro- 
file, so for instance a flat density proflle (7 = 0) would yield Ro 



Figure 1. Cartoon representing a star on a circular orbit in the 
combined potential of a black hole and a background galaxy. The 
black hole motion of amplitude -Ro runs parallel to the horizontal 
X-axis. The net force F acting on the star may be decomposed in 
a radial component F^ and an x-component. 



in the range 0.3 to 0.6 pc. 



changes sign when the phase ip shifts to ip + n. Thus when- 
ever the stellar phase-space density is well sampled and all 
values of if : [0, 2tt] are realised with equal probability, half 
the stars receive mechanical energy {SW > 0) and half give 
off energy {SW < 0). In other words, stars in the first quad- 
rant will exchange energy with those in the third quadrant 
of a Cartesian coordinate system. (Similarly for those in the 
second and fourth quadrants.) By construction, the black 
hole neither receives nor loses energy but merely acts as a 
catalyst for the redistribution of mechanical energy between 
the stars. Our goal, then, is to explore the consequences of 
this mechanism quantitatively for realistic stellar distribu- 
tion functions. 

We begin with an analysis of star-BH orbit coupling in 
a harmonic (uniform density) galactic potential (§2). While 
this choice may appear artificial and an over-simplification, 
it circumscribes all latitude allowed by uncertainties in the 
spatial distribution of stars within the black hole infiuence 
radius. Furthermore, the basic mechanics is more tractable 
for that case. This is then extended to the case of a loga- 
rithmic potential (§§3 and 4). To cover a wider range of pa- 
rameters, we explore with a response code the evolution of 
individual orbits in the time-dependent potential. We show 
that black hole motion shapes up the energy distribution 
function, as well as the line-of-sight velocity, which we mea- 
sure as rms deviations from expected values. The magnitude 
of these deviations rise monotonically with the amplitude of 
the black hole's orbit, and its mass. Finally, in §5 we dis- 
cuss some applications and explore possible extensions to 
our analysis. 
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2 Circular orbits in a harmonic potential 



We start with the case of a star initially on a circular path in 
a background harmonic potential. The star's orbit for that 
problem is obtained by solving the equations of a decoupled 
oscillator; these read in Cartesian coordinates 



(4) 



where 



iu = V47r Gp/3 

is the harmonic angular frequency in an axially-symmetric 
galaxy of uniform density p. Adding a fixed black hole of 
mass Mbh at the centre of the coordinates preserves the sym- 
metry of the force-field: the motion may still be described 
by Q but with a different angular frequency lu' > ui. Each 
circular orbit of the harmonic potential maps to a circular 
orbit in this new potential. The aim, then, is to find out what 
happens once the black hole is set in motion, so breaking the 
symmetry of the force field. 



W 



V$d« 



(7) 



The only non-vanishing contribution to the integral comes 
from the non-axially symmetric term of the force field; we 
find: 



W = - 



GMhh R 



GMbh R 



sin 9 rdO 



(8) 

For circular motion, we may always write rA6 = VcAt with 
both circular velocity Vc and radius r held constant. The 
ratio Vc/ ijjjr) satisfies 



\ru) / 



= 1 



^1 + 



Mg(< r) 



(9) 



where A4 £ [l,oo[ and Mg is the integrated galactic mass 
inside the orbit of the star. We have A4 = 2 when r — Ri,h by 
definition, and limr^oo = 1. To progress further requires 
a specific form for R{t). Solving (jljl for a radial orbit down 
the X-axis we find 



2.1 Co-planar, radial BH motion 



R{t) = Ro s\n{Lut + < 



(10) 



We consider a two-dimensional system so all orbits are copla- 
nar. Let the position of the black hole be denoted R and we 
take Mbh ^ The black hole obeys the same equations 
Q for any value of R not exceeding the uniform-density core 
of the model potential. (A Bahcall-Wolf cusp would soon 
develop around the centre once the black hole has settled 
there; the harmonic potential would remain largely unper- 
turbed until that happens.) We define the black hole radius 
of influence Rbh such that 



Mb, 



iTvpRtn/S 



(5) 



which we will use again later in a study of the logarithmic 
potential. As a flrst important case study we set the black 
hole on a purely radial path down the x-axis: R = Rx (from 
here onwards a hat denotes a unit vector). Equation Q 
gives an expression for the mechanical work for the case 
of a one-dimensional force component. We here develop a 
second approach based on a limited series expansion of the 
potential. The full potential 



simplifies to 



$(r,t) 



GMb. 



R' 



^+0([7?/rf; 



(6) 

when we take r > R and truncate a Taylor expansion to 
leave out all terms of order higher than quadrupolar. In this 
limit the force acting on the star is easily obtained by differ- 
encing ([6| with respect to r. A closed integral through one 
revolution on the circular stellar path I yields the net work 
done on the star by the time-dependent potential: 



where Ro , (po are the amplitude and phase of the black hole 
orbit. The initial conditions are completely specified if we 
pick the phase of the stellar orbit such that the azimuthal 
angle 6 = when t = 0. On integrating (jsl, we obtain 



W = - 



Ro GMg{< r)M 



27r 



(11) 

Written in this way we recover a direct proportionality be- 
tween W and the amplitude Ro (and hence W = for zero 
motion, as expected). W assumes both positive and neg- 
ative values according to the phase 4>o- There can be no 
correlation between cj)o and the star's orbit, and hence at a 
given radius r (or, equivalently, binding energy E) as many 
stars receive a positive energy increase as those that receive 
a negative contribution. This impoverishes the occupation 
level at r (or, E) and creates a hollow feature in the stellar 
distribution function. 

The work VK is a periodic function of (j>o ■ [0,27r]. In 
this interval, W vanishes whenever any of the following con- 
ditions is met: 



2-K \" TT 
TT — 2(j)o ) 2 



M = 



2-K 



Stt — 2(\>o y ' 2 ^ 

2-K y Stt 
57r-2<^J ''^°> 2 



(12) 



The full range of values of M are covered in the quadrants 
— 7r/2 < (j>o < 7r/2 and 7r/2 < 4>o < 37r/2. As a function of 
cf>o we find as many points outside as inside the radius of 
infiuence M = 2 (Fig. [2|. For an orbital configuration such 
that Ro < r <^ Rbh a large ensemble of solutions to ( |12[ | will 
cover a wide range of values of A4 > 2. This will be the case 
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M[] M[] 



Figure 2. Loft-hand panel: the net work W as function of the dimensionless parameter A4 defined in Note that W has been 
renormalised for better clarity. Right-hand panel: the angular momentum dL accrued during one revolution as function of Ai. The curves 
are labeled with the phase angle tpo- Note that the radius r would increase from left to right on the figure. The vertical straight line at 
Ai = 2 corresponds to r = Rth, while Ai = 1 when r — > oo. 



during the final stages of a black hole settling at the heart of 
a galaxy since then Ro 0. It is interesting to investigate 
whether configurations satisfying ( |12[ ) yield islands of stabil- 
ity or attractors in phase space. Fig. [2] shows the net work 
done on a circular stellar orbit by the black hole as func- 
tion of M for different values of (f/o- To determine whether 
an orbit will be trapped around a point for which W = 0, 
it is sufficient to look at the sign of W in relation to A4. 
Orbits with TVf < 2 lie outside the influence radius Rth at 
large distances. A gain in binding energy would shift the 
orbit to smaller radii requiring < 0. Similarly, an orbit 
for which M ^ 1 would shift to higher energy whenever 
W > 0, and so move outwards to larger distance (decreas- 
ing M). Thus curves of W for which dW/dA4 > when 
W = would trap orbits, otherwise not. We note that the 
migration induced through W impacts directly on the angu- 
lar momentum of the star. This is most clearly seen from the 
torque T = dL/dt = x r = -GMth/r^R{t)smez. On 
integrating over one (stellar) revolution it is easy to show 
that |dLj oc rW. At a given radius, the sign of dL is a func- 
tion of (f>o alone (see Fig[2jb]). For a well-mixed systems, half 
the stars will gain angular momentum and move outwards. 
Since the background galactic potential is taken to be axi- 
ally symmetric, the distribution of angular momentum at r 
can then be traced back directly to the quality of the black 
hole orbit (i.e., its amplitude and frequency). 

It is interesting to map out the minimum and maximum 
values of W achieved as a function of (/>o. An extremum of 
W occurs when the angle (po satisfies cos(27r/-v/!M + 4>o) ~ 
cos(0o) for any constant Ai (or equivalently, r). Fig. [s] 
graphs the results as a function of the phase angle (j>o- The 
results were obtained by setting a range of : [1, 100] for 



the mass variable. Recall that ^ oo as r ^ 0. The two 
curves on Fig. 3 coincide if we perform a reflection through 
W — together with a shift of 4>o by vr radians (180°). The 
figure shows that when a (positive) maximum is large, the 
corresponding (negative) minimum is small, and vice-versa. 
It follows that for a finite interval of the phase angle (po, for 
instance, all stars which fall in that interval may suffer a net 
positive (or, negative) intake of energy, irrespective of their 
orbital radius. Therefore for adequate initial conditions, the 
energy input from the black hole may leave a signature in 
the distribution function of the stars robust against any bias 
that might be attributed to domain decomposition (sam- 
pling by radius, energy, etc). By contrast, stars for which 
max(M^) « I min(W')| would suffer essentially no effect from 
the black hole. Initial conditions such that (j)o ~ n/3 radi- 
ans (« 52.7°) or 00 ~ 47r/3 (233°) are possible examples (cf. 
Fig.|3f. 

By virtue of ([9| and ( 11 1, we find in the limit of large A4 
that the work W oc M^^^ — > oo. Hence W diverges as r — > 0. 
This should not come as a surprise since in that limit the 
orbit is Keplerian around the moving black hole, and hence 
it can not remain close to circular about the centre of co- 
ordinates, as we have implied so far. Our development will 
therefore break down when the potential is dominated by 
the central point source. With M = 100, the largest value 
considered here, the model galaxy contributes 1% of the dy- 
namical mass only. For the MW galaxy, this would translate 
to a radius around Sgr A* of ~ 0.05 pc (Genzel et al. 2003). 
Stellar collisions are predicted to be important on a scale of 
0.1 pc (Yu & Tremaine 2003; O'Leary & Loeb 2006; Merritt 
2006). Therefore, if Milky Way data serve as a test case, we 
should always set A4 100. 
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Figure 3. This graphs the maximum (solid) and minimum (dash) 
value of the work W as a function of the phase angle (j>o (in 
radians). The two curves coincide when reflected through W = 
and shifted by Sij>o = i" rad. 



A final comment concerning tiie regime 3> 1. The 
definition (|9| is tlie ratio of tiie star's angular frequency to 
the galactic angular frequency oj. Thus whenever M ^ 1, 
the star revolves rapidly around the black hole. Such trapped 
orbits remain circular to a good approximation when viewed 
in the reference frame of the black hole. (The eccentricity 
e = is an adiabatic invariant.) Orbit trapping and reso- 
nant relaxation of non-circular Keplerian orbits have been 
discussed by various authors (e.g., Tremaine 1995; Ranch & 
Tremaine 1996; Zhao, Haehnelt & Rees 2002; Merritt 2006). 

2.2 Circular black hole motion 

2.2.1 Single black hole 

It is straightforward to extend the case of radial black hole 
motion to one where the black hole is on a circular orbit. The 
results of §2.1 are independent of the sense of rotation of the 
stellar orbit. The non-zero black hole angular momentum 
Lbh — R X p breaks this invariance. The work W will differ 
when stellar and black hole momenta are aligned {Lbh ■ > 
0) or anti-aligned {Lbh ■ L^ < 0), a shown by e.g. Toomre 
& Toomre (1972) in their classic study of interacting spiral 
galaxies. 

Consider a black hole on a clockwise two-dimensional 
circular orbit of radius Ro, 



R{t) = Ro {sin{ujt + (t>o)x + cos{ujt + (t>o)y] 



(13) 



in a self-evident extension to ( 10 1. Repeating the same steps 
that led to (|6| will give an extra term for the y-component 
in the potential which is otherwise identical. When comput- 
ing the work ^ we may now distinguish between anti- and 



clockwise stellar orbits with the notation d^ — ±Vc/rdt. 
The integration yields 



W 



Ro GMgM 



(iTl) sin[ 



27r 



+ < 



(14) 

This equation shows that aligned orbits double their energy 
intake from the black hole, while anti-aligned ones see a 
net cancellation. The black hole would, therefore, introduce 
anisotropy in an initially isotropic stellar velocity distribu- 
tion function. This is not unlike the bar amplification process 
proposed by Lynden-Bell (1979): the angular frequency uj set 
the rotation speed of a constant-magnitude quadrupole, as 
in a barred galaxy. However here the perturbation to the ax- 
isymmetric galactic potential is not the two-fold symmetric 
m — 2 mode, but the m — lopsided mode. 



2.2.2 Binary black hole 

The case of a binary black hole of constant separation and 



centered at the origin of coordinates is derived from ( 14 1 



through a thought experiment. We imagine that the binding 
energy of the binary is large and so neglect the background 
galactic potential. We consider the total contribution of the 
binary to W as a sum of two single BH's on circular or- 
bits. The phase (j>o of one black hole (say, the secondary) 
is shifted by tt radians with respect to that of the primary. 



On inspection of (14 1 it follows that whenever the product 



RoM is the same for each BH, the total work W must van- 
ish for any phase angle ipo. This will be the case only when 
the binary as components of equal mass. Since there is no 
reason to presume identical masses in general, we would find 



W 



Ro,2 — Ro,l GMg ,^ ^, 



■ r 27r 

sm , + I 



- sm 



where we have used Ro.i = Ro,2 Mbh,2/Mbh,i. The limit 
where Mbh,2 ^ Mth.i reduces to the case of a single BH 
(the secondary) in orbit in the axi-symmetric potential of 
the primary at rest at the origin of coordinates. Such a situ- 
ation might occur when a swarm of intermediate mass black 
holes revolve around a massive hole, presumably the result of 
repeated coalescence. O'Leary & Loeb (2006) have recently 
explored the scattering of stars in such a cluster of massive 
objects. 



3 Case study: the logarithmic potential 

The coupling between black hole motion and orbits in the 
harmonic potential is indicative of trends that may develop 
in more realistic potentials. Here we recast our problem in 
the framework of the logarithmic potential, which we write 



+ yVq' + Rl 
Rl 



(15) 



with Vo the constant circular velocity at large distances, and 
g ^ 1 is a (dimensionless) shape parameter. The radius Ro 
defines a volume inside which the density is nearly constant. 
Thus when r <^ Rc we have once more harmonic motion of 



6 Boily, Padmanahhan & Paiement 



angular frequency uj — Vo/Rc- If we let q 
u = r/Rc, the volume density p reads 



ATvGp{u) = V^$s = 
and the integrated mass Mg{u) 



2^ -"o 3 + 



Rl (1 + V? 



1 and define 



(16) 



MJ< u) = 



V^oRc 



G u2 + 1 



(17) 



The mass Mg{u 2> 1) oc u diverges at large distances, how- 
ever this is not a serious flaw since we will consider only the 
region where u ~ 1. The mass Mg{u = 1) = vlRc/2G fixes 
a scale against which to compare the black hole mass Mi,h- 
Since the black hole will orbit within the harmonic core, we 
set 



M, 



'bh = rubh- 



2G 



(18) 



with < rhbh ^ 1, and 



M{u) = 1 + rhbh 



1 + 

7,3 



bears the same meaning as before. The core radius offers a 
reference length to the problem. The position and velocity 



of the black hole at any time follow from (10 1, where the 



amplitude is set by fixing the dimensionless number Uo — 
Ro/Rc and the angular frequency uj is 

_ ^ 
Rc 

Our goal is to quantify the time-evolution of a large num- 
ber of orbits in the combined logarithmic and black hole 
potentials. If we pick parameters such that 

m* < A4bh < M5(max{u}) 

then we may neglect the collective feedback of the stars on 
the black hole and galactic potential and study only the 
response of individual orbits evolving in the time-dependent 
total potential. This approach will remain valid so long as 
the response of the stars are relatively modest. The time- 
evolution of orbits was done numerically using a standard 
integration scheme, which we describe below. 



3.1 Equations and numerics 

The energy Ej of a star is computed from 



Ej = ^v^ + ^g{r) + <i'bhir,t) 



where 



^bh{r,t) = 



GMbh 
\r-R\\ 



(19) 



(20) 



is an explicit function of time through (10 1; from here on- 



ward we will write {xbh, yth) for the coordinates at time t of 
the black hole in the a;y-plane. The time-derivative 



Ej = dt^bhir,t) j^O 

and hence energy is not a conserved quantity. We use this 
fact to define a set of six first-order differential equations 



dt 



(r,Ej,v,Ej) = {v,Ej,-V<i>,Fj) 



(21) 



where the time-derivatives are computed in the usual way, 
and we find for a black hole orbit confined to the harmonic 



Fj 



^ \\r-^R\\5 ([^^''^ ~ ^] [xT^x] + [ybh - y] [yTh^y]^ 



GMb, 



-R\\' 

GMb, 



(xlh + ylh - ixbh 



yyth 



([xbh ~ x\xbh + [ybh - y]ybh) 



(22) 



\\r-R\\3 

The expression for Fj admits a simplification when the black 



hole is set on a radial orbit (10 1 but note that (22 1 is not 



invariant to a swap of x for y and vice versa when ybh ~ 0. 



3.1.1 Compact kernel 



The potential ( 20 1 is singular when r — R which introduces 



large errors in the integration. To alleviate this we redefine 
( 20 1 using a compact kernel, effectively smoothing over the 
singularity. Let Sbh be a constant length and write 



^bh(Tl,t) 



-^$(7^,t), 

£bh 



(23) 



whenever TZ ^ 1 where TZ — | |r — i?| |/e(,fi. We pick a compact 
kernel which minimises force errors at 7?. = 1 (Dehnen 2003) 
and define 



$(7^) = 1 



i(l-7^^) + |(l 



(24) 



This last equation fails when TZ > 1, however this is of 
no concern since the gradient is continuous at 7?. = 1 and 



matches exactly the one derived from ( 20 1 at that radius 



Integration of equations (211 with respectively (20 1, or (23 1 
and (1241), when 7?. > 1 or ^ 1, poses no particular difficulty. 



though Eq. 22 takes another form inside TZ < 1 (see below). 



3.1.2 Choice of units & integrator 

For convenience we have chosen scales for the background 
potential such that G — Vo = M = 1. Borrowing from the 
case of the MW black hole, we set a kernel length Sbh — 2 x 
10~^ which will wipe out all high-deflection angle collisions, 
i.e., those due to orbits with little angular momentum. 

We have used an explicit fourth-order time-adaptative 
Bulirsch-Stoer integrator taken from Press et al. (2002) for 
solving 1 21 1. We have performed a series of tests with a 
sta tic p otential by setting e.g. rhbh ~ 1 and Uo — Ro = 
in (10 1. With these parameters Ej = 0, and we checked 
that a precision of 1 : 10^'* is maintained for a runtime of 
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200 units. In particular we validated Eq. [24] by integrating 
radial stellar orbits running through the black hole, both 
along the x- and y-axis. We also integrated a circular orbit 
at the edge of the kernel, ov TZ= 1, and found no indication 
of a drift in energy or any kind of random fluctuations. 

The situation is less glorious when integrating with fi- 
nite black hole motion. To see why, let us write down Fj in 



(211 with (10 1 and TZ < 1. Some straight forward algebra 
yields 



Fj 



3—5— (x- 



3^2 

2^ 



At 



(x- 



XbhjXbh ■ 



(25) 



Since the relative distance between the star and the black 
hole is < e, we find that Fj < 0(l/e^) which gives O(IO^) 
for the parameters chosen. As such this would not be prob- 
lematic, however the time-steps required to maintain an ac- 
curacy of 1 : 10^* for a typical integration time become tiny, 
and the computer run-time, prohibitive. We found a practi- 
cal solution to this problem, by imposing that the quantity 
e^Fj be integrated to a precision of 2 : 10^^, so that Fj is 
known to six significant digits (e ~ 10"'^). Whenever this 
condition was not met, we only included the orbit in the 
analysis up to that point in time, after which it was ignored. 
This situation occurred relatively seldom, and affected some 
3 — 4% of cases at most. This was so, for example, when the 
initial configuration either started out with many stars on 
near-radial orbits and small-amplitude black hole motion; 
or, when the black hole was allowed to flirt with a large 
number of stars by covering a distance comparable to or 
exceeding its radius of influence. 



3.2 Initial conditions 

We limit our exploration to the case of coplanar motion. 



The volume density ( 16 1 stretches to infinity and yields a di- 



vergent integrated mass. Since we are only interested in the 
central-most volume, we decided to keep only stars that re- 
main inside a given radius. (A selection by energy Ej would 
be equivalent.) Our model calculation of §1 had shown that 
the response of the star to black hole motion is a strong 
function of the ratio of their orbital periods. The orbital 
period of the star is oc 1/ \fGp, where p is the mean volume 
density inside the semi-major axis of the orbit. These two 
observations combined suggests that we only include orbits 
out to where the density varies most rapidly. We chose 
a truncation radius such that ^ e, close to the value 
2.80(7) at which din /o/d In r = —2.101.. reaches a minimum. 

We now specialise to the case of circular orbits in the 
axisymmetric potential $5 by setting q= 1 in ( 15 1. Positions 



are attributed by Monte Carlo method using the density as 
probability distribution. The square circular velocity at each 
radius is 



tt2 + 1 



"2ir 



(26) 



tal angular momentum of the stars is zero to within root-n 
noise. The energy may be written as 



— — = \u\u + 1 



+ 1 



TTlbh 

2u 



(27) 



This expression is easily differentiated to yield an analytical 
form for E'j = AEj jAu which is the density of states of stars 
of energy Ej at u. Since we are only interested in coplanar 
orbits, the mass element drawn from (|16|) is 



<5M = 2-KRip(u)uAu = f{u) du = f{u) 



dEj 



f{Ej)dEj 



where the energy distribution function f{Ej) is known in 
parametric form. 



G mbh(«2-f l)2+4«3(7i2-(-2) 



(28) 



This equation shows that when rhbh = (no black hole) we 
find f{Ej) — > constant in the limit u ^ 0; and f{Ej) oc u'^ 
in the same limit when rhbh 7^ 0. Thus the bulk of the stars 



avoid the central black hole. Equation ( 28 1 will be helpful 



when assessing the noise level of the response of the stars to 
black hole motion. 



4 Results 

We have until now fixed the gravitational constant G = 1 
and velocity scale of the galactic potential Vo = 1. The initial 
conditions require further that we fix the black hole mass 



parameter rhbh in ( 18 1, black hole's amplitude of motion Uo, 



and the core radius. Re- The total mass (17 1 will integrate 
to 1 up to Ut — e if we set R^ ~ 0.46. This fixes all scales 
in the problem, and we note that « half the stars lie inside 
Rc, and an equal fraction outside. The black hole's radius of 



infiuence is obtained in terms of rhbh from equating ( 17l to 
( 18 1. The result is shown on Fig. [4] We use this relation to 
set a more stringent constraint on the motion of the black 
hole by imposing that it orbits always within its radius of 
influence, an improvement on our initial ansatz that Uo < 1, 



since ( 10 1 is a better solution to the black hole's orbit closer 



and the sense of motion chosen randomly so that the to- 



to the origin. 



4.1 Reference case, rribh = 0.3 

We set up a reference case inspired by MW data which will 
guide us through our exploration of parameter space. Since 
the MW black hole lies close to the galactic centre at a 
velocity ~ 2 — 4 km/s this suggests that we focus on cases 
where the black hole remains well inside the central core. We 
pick a black hole mass equal to 30% of the core mass, rhbh ~ 
0.3, and set an upper limit of Uo = 0.66 on its amplitude of 
motion, when it would exceeds slightly its influence radius 
of ~ 0.57 (Fig.|4|. 

To see how orbits respond as Uo is increased from zero, 
we draw on flgure [5] Poincare sections of a single orbit 
for four values of Uo : 0,0.15,0.33 and 0.66. The case of 
a flxed black hole is shown on Fig. [5|a), when the star 
describes a circular orbit of radius r = 1.78 which is > 3 
times the influence radius of the black hole. The orbital 
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Figure 5. Poincare section {x,x) at y = (left panels), orbit (middle) and radii and energy (right) for a single stellar orbit in the 
logarithmic potential to which we added an rhbh = 0.3 bla<;k hole. The panels to the right display the binding energy of the star (solid 
line) along with the distance r* of the star to the black hole (dash) . Modulations in energy match one to one variations in r* . The black 
hole orbit rf,h (dots) is also displayed for comparison. The four rows show the orbit for different values of Uo'- a) 0.0, b) 0.15, c) 0.33 and 
d) 0.66. 
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mass parameter m-tilde [ ] 

Figure 4. Influence radius in units of Rc as a function of the 
mass parameter mi,h defined in llSjl. In the limit rrt;,;, —> the 

1 /3 — 
radius oc rises rapidly. In the neighbourhood of mj^ = 1, 

the radius oc mj^; the straight line is the curve mf,^/2 + |. 



Table 1. Initial conditions of the numerical orbit integrations. 
The influence radius ri^f^ deflned in |5]l is given in computational 
units. The potential l |15| is deflned in units such that G = Vo = 1. 
The core radius Rc = 0.459, and the truncation radius rt ~ 2.53 
gives an integrated mass ( |17| = 1. 

Circular stellar orbits ('cold') 



Name 


rhbh 


Uo 


rbh 


Comment 






[Ro/Rc] 






CI 


0.30 


0.00 


0.258 


Static 


C2 


0.30 


0.07 


0.258 




C3 


0.30 


0.15 


0.258 


Reference case 


C4 


0.30 


0.21 


0.258 




C5 


0.30 


0.30 


0.258 




C2s 


0.15 


0.15 


0.197 


Shadows C2 


C3s 


0.15 


0.30 


0.197 


Shadows C3 






'Warm' or 


'Hot' runs 


Name 


rhbh 


Uo 


rth 


Comment 






Ro/Rc 






Wl 


0.30 


0.15 


0.258 




Wlc 


0.30 


0.15 


0.258 


circular BH orbit 


HI 


0.30 


0.15 


0.258 




W2 


0.30 


0.21 


0.258 





period — 2nr/vc — 11.18 and the integration was for a 
total of 200 time units (18 revolutions). The middle- and 
right panels show the orbit and the star's energy (solid 
curve) and distance to the black hole (labeled r^,, dashed 
curve), respectively. Figjsjc) illustrates the situation when 
the black hole motion has amplitude Uo — 0.33 ~ 0.58 x its 
influence radius. The star's orbital radius now varies from 



a minimum of approximately 1.65 and up to 1.85, a gap of 
~ 10% compared to the circular orbit; the same conclusion 
applies to the cycles seen in binding energy (right panel. 
Fig. [sjc]). We note that the modulations in E match one 
to one the profile of r*, which is an indication that a 
strong coupling is still operative even for stars orbiting well 
beyond the black hole's radius of influence. The scatter 
seen in the Poincare sections of both Fig. [sjc) and (d) 
confirms this view, and suggests that the original circular 
orbit becomes mildly chaotic (no loop or resonant structure). 

The response of individual orbits to black hole motion 
would leave a measurable trace only if their signal rises above 
the background noise of other orbits that may be otherwise 
affected. We sought out a relation between an unperturbed 
and time-independent distribution function and the noise 
level of a discrete realisation of that function with N bodies. 
The reference distribution function is given by \28\ which we 
discretised using 100 equal-size bins of width AE ~ 0.07. 
This curve is plotted as a histogram on Fig. |6] We then 
measured the root-mean square differences with the ana- 
lytic distribution function by drawing different numbers of 
orbits in the range 10,000 to 100,000. The results are shown 
for four values of TV on Fig|6];a-d). We find the rms differ- 
ences drop to ~ 4 — 5% already for TV > 40, 000, comparable 
to Poisson noise (roughly l/\/400 ~ 5% fluctuations per 
bin). For completeness, we also plot the relative differences 
in percentage at each bin of energy on the rectangular frames 
below each panel. This relative energy error is dominated by 
low- number statistics and becomes very large when i? < — 1, 
which, for our choice of parameters, corresponds to a radius 
u ~ 0.074 enclosing ~ 0.02% of the total mass (expected 
Poisson noise of ~ 20%). The error made in dropping orbits 
below that level of energy from our analysis is of the same 
order. Furthermore, stars that are that close or closer to the 
origin would be trapped by the black hole gravity and re- 
main on high-velocity Keplerian orbits around it. This would 
remain true even if the black hole were set in motion at a 
comparatively small velocity. Such orbits are not the focus 
of this work. 

Figure [7| compares the energy distribution function of 
a set of 93617 orbits when the black hole is set in motion 
with amplitude Uo = 0.33, to the initial distribution func- 
tion, when the black hole sits at the origin. Recall that this 
latter distribution would be time-independent. The orbits 
were all integrated for t = 20 unit of times, which corre- 
sponds roughly to 10 revolutions at the edge of the sys- 
tem. The energy axis is once more split in 100 bins, and 
we find once more an rms noise level of ~ 4%. Let us take 
as one standard deviation a difference of 5% with the ana- 
lytic function ( |28[ ). The signature of black hole motion seen 
on Fig. [7|b) leaves four peaks of more than three stan- 
dard deviations, and ten with one or more standard devi- 
ations, all to the left of _E = — 1. The most significant peak 
at £ ~ 1.31 has an amplitude of -1-29%, or five standard 
deviations. This corresponds to a circular orbit at radius 
u — 1.169 (r — uRc — 0.5355 units), which is ~ 2 times 
the black hole radius of influence (cf. Fig.|4| and 3.57 times 
its amplitude of motion, Uo- This is clearly an indication of 
energy exchanges through beat frequencies, i.e., resonances. 
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Figure 6. The energy distribution function (d.f.) for four discrete realisations with (from [a] to [d]): A'^ = 10,000,20,000,40,000 and 
80, 000 orbits. In each case the span in energy was divided in 200 bins of equal size. The analytic curve | |28| l is displayed as the more 
regular histogram on each panel, and the root-mean square differences given. The rectangular frames gives the relative differences (in %) 
for each bin. 



4.2 Orbital resonances 

To gauge the importance of orbital resonances we identify 
first the radius and energy of commensurate orbital periods. 
If we call uj' the orbital angular frequency of a star on a 
given orbit of energy E, then we need solve the equation 



1 rribh 
+ 1 2u3 



(29) 



where v = uj/ui' as in §1, which requires solution for all 
prime integer ratios m/n. The above equation could be set 
in terms of A4 defined in ^ as done in ( 12 I in the case of 
a harmonic potential. Instead we solve for E,u from (291, 



and classify the result as a Keplerian resonance when the 
corresponding value for M > 2, which will always be the 
case when m > n, and a Galactic resonance for all other 
cases m ^ n. Table [2] lists the results for a broad range 
of values of m : n. Remarkably, the energy level of Galac- 
tic resonances match almost exactly the apparent nodes in 
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Figure 7. Energy distribution for two configurations: A) static black hole at the centre of coordinates; B) black hole on a radial orbit of 
amplitude Uo = 0.33, after 20 time units of evolution. In both cases the mass parameter m;,;, = 0.3. The energy d.f. averaged over five 
output times is shown on the right-hand side for comparison. The smooth histogram on all three panels was constructed using l|28|. 



Table 2. Orbital resonances defined in \29\ for several values of 
the commensurate ratio m : n. Results for m S^n (or, M 1) are 
labeled Galactic, otherwise they are labeled Keplerian [Ai < 1). 



Galactic 




u r 


2i?/^2 


m : n 








1;1 





748 


0.602 


2:3 


1 


292 


1.490 


3:5 


1 


479 


1.744 


1:2 


1 


849 


2.178 


2:5 


2 


390 


2.691 


1:3 


2 


918 


3.097 


Keplerian 




u r 


2i?/^2 


m : n 








11:10 





663 


0.443 


6:5 





598 


0.317 


4:3 





532 


0.187 


13:9 





489 


0.101 


3:2 





471 


0.064 


5:3 





427 


-0.031 


2:1 





364 


-0.171 


5:2 





304 


-0.320 


3:1 





265 


-0.433 



the energy distribution function seen on Fig. [7|b). This is 
in close agreement with our analysis of §1, when we argued 
that commensurate values of u would give no net work after 
an integer number of black hole revolutions. This is not the 
whole picture, however, since the binding energy of individ- 
ual orbits is constantly changing in time, with as many stars 
gaining energy as those losing energy. Hence one may think 
of the energy d.f. of Fig.jTj^b) as a standing wave modulated 
by small sinusoidal modes propagating at a finite pattern 
speed. A clearer picture emerges once we compare Fig. [7|b) 
to a time-averaged d.f. for the same system. On the right- 
hand panel of Fig. [7] we graph the average of five energy 
distributions functions sampled over a time span of t = 13 
to 17. This corresponds to A/uj — ARc/vo ~ 2 full black hole 
oscillations. In total 468115 orbits were put to contribution. 
Comparing this curve to the one displayed on Fig. [7|b), 
we find fewer peaks exceeding one standard-deviation. The 
smoother appearance of the distribution function supports 



the interpretation of sine-like oscillations on Fig. [7|^b) as 
transitory features. Thus we expect phase-mixing to erase 
such features on a dynamical time-scale. By contrast, the 
two broad peaks at 2E/v^ « 1 and —0.5 remain, their am- 
plitude hardly dented by the time-averaging. These peaks 
should, therefore, leave observable features in kinematic- and 
density maps. 

This intuition is confirmed, at least partially, by a 
Poincare section of {x,Vx) in the plane y = 0. On Fig. [s] 
we graph the surface of section of 500 orbits, each evolved 
for twenty time units. The number of points varies between 
orbits from ~ 10 and up to 200, according to the period. We 
observe large but localised scatter in the velocities as a func- 
tion oiux = x/Rc- The vertical full line at li^, ~ 0.3 indicates 
the amplitude of motion of the black hole, while the dash 
lines are the location of resonances listed in Table [2] The 
scatter decreases rapidly as we move to large radii. Specifi- 
cally there is a sharp drop as we reach beyond the 1:1 reso- 
nance, and thereafter significant scatter is centered around 
Ux ~ 1.6 and 2.5. These values of Ux correspond to energies 
of 2E/vo ~ 1.9 and 2.8, respectively, matching the features 
seen on the energy distribution function (cf. Fig. [7|. The 
broad peaks seen on graphs of the energy d.f. falls inside the 
1:1 resonance and are lost in the scatter on Fig. |8] 

4.3 Surface density, velocity maps 

Maps of the surface density and velocity field are of inter- 
est. The configurations are isotropic initially when the black 
hole starts off at the centre of coordinates. At later stages 
neither the density nor the velocity fields respect this initial 
property. We toyed with the idea of plotting both surface 
density E and velocity field in cylindrical coordinates cen- 
tered either on the system's centre of mass, or the black 
hole. This turned out to be useful only when both coincides, 
which will only occur when the black hole is on a radial orbit. 
Instead, we opted to map out both quantities on a uniform 
Cartesian grid. This has the advantage of an unbiased lin- 
ear resolution over all space and is identically suited to any 
type of black hole orbit (radial, circular or other). The sur- 
face density is obtained at any time by a simple count-in-cell 
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Figure 8. Surface of section in the j; = plane showing in 
units of Vo as a function of = x/Rc- The dashed vertical lines 
indicate resonances listed in Tabled 



(CIC) technique. No smoothing or averaging of neighbouring 
cell has been performed. Density profiles and their signifi- 
cance should be interpreted with due consideration to root-n 
noise: a typical grid would have 30 x 30 mesh points, and 
hence a mean count per cell of at least > 4 10^/900 ~ 44, 
which translates to relative fluctuations of 15%. In practice 
we have used on order 10^ orbits, and so the noise level al- 
ways falls in the range 10% - 15%. As we will see, the density 
fluctuations that we measured at times exceeded 60% of the 
reference initial profile, giving a signal-to-noise ratio of at 
least 4 and perhaps as high as 7. 



4-3.1 Flat-fielding the velocity map 



The initial velocity 1 26 1 is known at any point in space, how- 



ever it is set in the centre of mass reference frame. Since we 
wish to map out velocities on a Cartesian grid, we must be 
cautious to compute the expected mean velocity in any cell 
for comparison. This poses a problem around the origin of 
coordinates, when the radius ~ a few grid cells only. Calling 



du the grid size, we find from (26 1 an absolute error on Vc 
of 



\SVc 



2u 



+ 



rrihh 



2 max{tt, e} 



(30) 



which becomes large when u < 5u. We then compute ( 30 \ for 



each orbit falling inside a given mesh and take the average 
square difference with the local circular velocity: 



n ^ ( 1 



+ \ 5Vc\/Vc 



(31) 



where the sum is over all n orbits inside the mesh at time 
t obtained by CIC. This gives a direct measure of the local 
dispersion as a result of black hole motion, and a reference 
map to eliminate noise when the black hole is fixed. For 



that case, we find using (31 1 residual errors not larger than 
1 : 10* or 0.01 %. 



4-3.2 Ttme sequence 

We graph on Fig.[9]the time-sequence of the surface density 
and dispersion a for 93607 orbits integrated over 20 time 
units. The figure shows two sets of two rows, regrouped to 
help match features seen in the density, to those seen in 
de velocity field. Initially the orbits are isotropic and circu- 
lar, which explains the two feature-less frames at the top- 
left corner of Fig. |9] At subsequent times, the plots show 
very pronounced and fast-evolving features, both in maps 
of the density and velocity. The scale of surface density was 
chosen to identify peak density enhancements of 60% when 
compared to the initial profile. The core radius Rc ~ 0.45, 
nearly twice the black hole influence radius of ~ 0.56-Rc, 
is displayed as the dashed circle. Inside and up to that ra- 
dius, the mass profile shows arcs, bubbles and other tran- 
sient features, all suggestive of unsteady, perhaps chaotic, 
orbital motion. Outside that radius, we find more steady, 
ring-shaped features which match the position of resonances 
(Table [|). 

By comparison, the relative velocity dispersion peaks at 
~ 25% of the local circular velocity Vc- Not surprisingly, the 
largest deviations in velocity nearly always coincide with the 
position of the black hole (red dots on Fig.[9|. The most re- 
markable features on these frames are the large dispersions 
measured well outside the core radius. A particularly strik- 
ing sequence runs from t = 15 to 17, when a large arc seems 
to close up on itself in the region x « — 0.75,j/ ~ 0. This 
large dispersion is found at a distance some 5 times larger 
than the amplitude of the black hole's orbit. If we scaled 
these features to the MW, then an anomalous local velocity 
dispersion would be expected up to from 3 to 4 pc away 
from Sgr A*. In conclusion, the black hole has a very strong 
impact on circular orbits, both in terms of spatial features 
and kinematics, up to ~ 2 to 3 times its radius of influence 
(see also Figlsl. 



4.4 Projected velocity field 

The line-of-sight velocity offers a direct way to measure the 
effect of black hole motion. Our strategy, then, is to apply 
a slit of some 100 mesh running across the x- and j/— axis 
(viewing angles of 0° and 90° to the black hole orbit, respec- 
tively), and compute both the integrated one-dimensional 
velocity and the square velocity dispersion at each mesh 
point. The expected velocity in projection vid{x or j/) = 
for isotropic initial conditions. The analytic calculation of 



the dispersion cy\d proceeds straightforwardly from ( 16 1 and 
(26 1. Inspection of several time frames, using either x or y 
as the line-of-sight, revealed no significant deviation from 
the profile derived for circular motion, at any point of evo- 
lution. Specifically, we looked for asymmetries in the dis- 
persion when comparing the two viewing angles, but found 
none (to within ~ 2x the noise level). Consequently we dis- 
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Figure 9. Time-sequence showing the surface density S and flat-fielded velocity dispersion a for a system with an m;,/! = 0.3 mass 
black hole on a radial orbit of amplitude Uo = 0.33 (Ro ^ 0.15). The frames show maps of S (first and third rows) and dispersion a 

(second and fourth rows) at times t = 0, 13, 14, 15, 16 and 17 units. The black hole is shown as a red dot, while the cross is the origin of 
coordinates. The black arrows are the local net angular momentum, which remains at the root-n level everywhere. 
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Figure 10. Time-sequence sliowing tlie line-of-siglit root-mean-square velocity dispersion, a, for the reference case of rhj,/! = 0.3 ~ Uo- 
A total of 100 bins were used for two viewing angles, down tlie x- and y-axis, in eacli case. At t = (panel to the extreme left) the 
analytic expectation with a black hole at the centre of coordinates is well recovered from 69773 circular orbits. The dashed line shows the 
theoretical expectations which drops to zero at the origin (all motion is orthogonal to the line-of-sight). The solid line gives the dispersion 
recovered from summing over all orbits. When t > 0, a differs from the expected value up to u^^y ~ Uo, the black hole amplitude of 
motion. At large distances, the broken curve suggests Poisson noise from small-N statistics. The thin dotted lines are the differences 
about the coordinate centre <j{u) — a-(—u) for each line-of-sight. 




Figure 11. Time-sequence showing the line-of-sight mean velocity for the reference calculation with black hole mass rrtjfi = 0.3 and 
amplitude Uo = 0.33. The velocity was measured when looking down the x-axis (solid line) and the y-axis (dashed line). The x-axis data 
show cyclical variations with a peak signal-to-noise ration inside u = 1 (or, r = ifc). 



carded runs of aid from analysis. (In §5 we comment on the 
dispersion in the context of fragmentation modes.) 

The situation with vid is more profitable. Fig. |ll| graphs 
the hne-of-sight velocity for three time frames of the refer- 
ence C3 calculation with rhbh = 0.3 and amplitude Uo = 
0.33. The results are shown for two viewing angles in each 
case at times t = 0, 14, 15 and 16. A total of 83697 orbits 
were binned on a mesh of 100 points. The Poisson fluctua- 
tions ~ 3% on average but are less than 1% inside |u| < 1, 
and ~ 10% near the end-points due to sampling effects. 
Looking down the y-axis orthogonally to the black hole or- 
bit, we find vid = everywhere to root-n noise, a result 
which confirms the intuition that the y-component of the ve- 
locity field preserves the initial symmetry through the z — x 
plane (Fig. 11 dashed lines). A sample of these three and 
two more frames gave very similar results for the root mean- 
square scatter \/ < 
mesh. 



0.002 when averaged over the 



When we switch lines-of-sight, we find large time- 



dependent oscillations of vid of a half-period ~ 2, close 
to half the orbital period — 2-kRc/vo — 2.9 of the black 
hole (Fig. 11 solid line). This highlights a close relation 
between the stars' angular momentum = xvy — yv^ 
and the phase of the black hole orbit: the torque F^ ~ 
GMbh/r'^ (i/*/r)(a;* - Xbh)z (where r = ||r* - r^hH is the 
distance between the star and the black hole) will be pos- 
itive for half the stars inside a given mesh centered on j/*. 
This torque, then, will give a boost to the momentum of 
stars when ■ Lz > (otherwise the torque will oppose 
the star's momentum, and decrease its magnitude). As a 
result the line-of-sight velocity does not integrate to zero 
and shows cyclic variations with Uy . Inspections of the same 
five frames as before using the a;-axis as line-of-sight showed 
that \vid\ varied from 0.008 to 0.034 with a root mean square 
scatter y^<Svfa> « 0.03. The run from t = 14 to 16 is an 
example of a sequence during which |«id|, averaged over all 
bins, goes from 0.032, to a minimum 0.004, and then back 
to 0.034. We stress that the trend with time is averaged over 
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Table 3. Results for the models listed in Table 1. 



Circular stellar orbits ('cold') 



Name 


rms {25E/vof 


rms (<5f 


)/fo 


Comment 




±4[%] 


X-axis 


y-axis 




CI 





0.015 


0.015 


t = data 


C2 


12 ±2 


0.022 


0.020 




C3 


29.2 


0.071 ±0.01 


0.042 


Ref. case 


C4 


36.9 


0.044 


0.026 




C5 


51 


0.038 


0.024 




C2s 


24 


0.025 


0.021 


Shadows C2 


C3s 


50 


0.041 


0.022 


Shadows C3 




'Warm' or 'Hot' runs 




Name 


rms {25E/vo)^ 


rms (<5f^^ 




Comment 




±4[%] 


X-axis 


y-axis 




Wl 


18.8 ± 1.6 


0.050 


0.015 




Wlc 


22.0 


0.041 


0.031 


circ. BH orbit 


HI 


28.5 ±3.5 


0.034 


0.014 




W2 


32 ± 1 


0.033 


0.016 





all bins; it is therefore of a much larger amplitude than the 
scatter seen inside \uy \ = 1. Alternatively, we may compute 
the ratio of \vid\ measured down one axis, to the values ob- 
tained for the initial configuration &tt = 0. We find a scatter 
^/<^vf^ ~ 0.015 averaged over all bins, but only ~ 0.005 
inside |w| = 1 owing to better statistics. If we refer to these 
data as noise, then the signal-to-noise ratio when goes 
through a maximum reaches 10 near u ~ at these times. 
The rms scatter of x-axis data is systematically ~ 3 x higher 
that of y-axis data, at all times (Table [3|. 

The impact of black hole motion on the distribution 
of angular momentum L appears clearly on a graph of the 
momenta distribution function. The panel to the top left on 
Fig. 12 (labeled 'cold') shows the d.f. at time t = 17 for 
the reference calculation, compared to the initial profiling 
(thin dashed line). The angular momenta are distributed in 
a highly symmetrical fashion about the L — axis. Each 
feature marking a departure from the initial distribution is 
matched pair- wise for the same value of \L\. Thus the sum 
of all stellar momenta remains constant. Fig. 12 (top, left- 



hand panel) illustrates the role of the black hole acting as 
a catalyst to transfer angular momentum, as well as energy, 
to the stars. 



4.5 Exploring different configurations 

4-5.1 Changing the velocity field 

Our choice of circular orbits has some bearing on the out- 
come of the calculations. In this section we assess to what 
extent features seen on Fig. |9] are specific to our choice of 
initial conditions. 

To that end, we first took the same set of orbits but 
modified initial velocities by ±10% in magnitude, so that 
the orbits were no longer circular initially. This exercise pro- 
duced similar features in the maps of density and velocity 
as for the reference setup, a hint that the initial response 
of the stars, and the general features, are not sensitive to 
imposing strict circular motion to the initial conditions. We 
then toyed with the idea of computing a velocity field consis- 



tent with the logarithmic potential based on moments of the 
Boltzmann equations (e.g., Binney & Tremaine 1987). Do- 
ing so however would have meant introducing a new energy 
d.f. and made comparisons with the case of circular orbits 
more difficult. Instead, we opted to keep the same energy d.f. 
and to change the velocity field by reorienting the velocity 
vectors randomly inside some angle chosen in the interval 
±9; here 9 — gives the original d.f. with all orbits circular. 
Two new configurations were setup, one with 6 = n/'i (giv- 
ing a cone of opening angle 45°), which we label warm, and 
a second with 9 = tt (fully random in azimuth), which we 
label hot. The distributions were otherwise unchanged from 
the reference model C3. The full list of model parameters is 
given m Table [l] 

The results are summed up graphically on Fig. |12[ The 
top three rows show the angular momentum d.f. (left-hand 
panels), the surface density (middle) and vid for two projec- 
tion axis (right-hand panels) in turn for the 'cold' configura- 
tion of circular orbits, the 'warm' and 'hot' initial conditions. 
All data shown on the figure were taken at time t = 17 so 
the black hole assumes the same position and velocity in 
each case. The rings seen in the surface density map of the 
cold run outside r = i?c ~ 0.5 (shown as a dashed circle) 
have disappeared, and only the strongest feature inside Rc 
remains visible for warm and hot initial conditions. Thus 
the filamentary structures seen on Fig. |9] are attributable to 
the strong response from circular and near-circular orbits. 
These would make up a small fraction of warm and hot dis- 
tributions. In these two cases, the number of low-L orbits 
is higher and more stars visit the central region on eccen- 
tric orbits. These scatter off the black hole effectively and 
acquire large angular momenta, a feature which can be mea- 
sured up from the swelling of the momentum d.f. for large 



12 



\L\ and the depletion of the d.f. around L = (Fig, 
left-hand panels, 2nd and 3rd row). The black hole motion 
has Uo — 0.15 translating to r = 0.069 and a radius of in- 
fiuence ~ 0.23. Therefore all stars with an initial angular 
momentum lower than r x i; « 0.34 would come within a 
distance < rf,h of the black hole in one revolutiorj^ Fig. 12 
shows this estimate in good agreement with the numerical 
computations. A fraction of stars concerned is about 12% 
in the 'warm' calculation, but close to 25% in the 'hot' run, 
when circular orbits have all but been wiped out. 

Although details of the density and angular momentum 
profiles are much affected in the new configurations, com- 
pared to the reference cold one, the same does not apply 
to the line-of-sight velocity. The response of stars for both 
warm and hot computations show a distinct signature of 
black hole motion in the sense that once again vid fluctu- 
ates significantly more when measured down the axis parallel 
to the black hole's orbit. The averaged scatter in both the 
cases compares well with the data for the cold configuration 
of circular motion, at least in the interval from it = to 
jit I = 1, or roughly two times the black hole radius of in- 
fiuence (Fig. |12| and Table [3|. In conclusion, the black hole 
motion still imprints the kinematics of both warm and hot 
configurations, albeit to a lesser degree than for the case 
when stars are on circular orbits. 



We have set r = 0.23 + 0.069 or u = 0.65 and v{u) from l26b 
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Figure 12. This graphs from left to right: the angular momentum distribution, the surface density, and the line-of-sight velocity for four 
models. The top row gives the results for the reference C3 calculation; the second row is the 'warm' Wl calculation; and the third row is 
the case HI of a 'hot' distribution (see Table[l]and §4.5 for details). All quantities were analysed after t = 17 time units of integration; a 
minimum of 93 000 orbits were used to sample each of the parameters. The last row at the bottom gives the solution for the case Wlc 
of a black hole set on a circular orbit. Note the strong m = mode on the map of the surface density. 
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2E/Vo2 




Figure 13. This graphs the energy distribution for four values of the black hole's amplitude of motion Uo = Ro/Rc- The dimensionless 
mass parameter m^^ = 0.3 in all the cases. On panel (a) we plot two solutions with = 0.07 (solid line) and 0.15 (dashed line); on (b) 
we set Uo = 0.21 (solid) and 0.30 (dashed). The smooth curve is the analytic solution (28). The dots shows a discrete realisation with 
69977 orbits and Uo = (axially symmetric potential). 



4-5.2 Changing the black hole parameters: scaling 



Equations (111, (17 1 and (181 may be combined to give a 



proportionality relation between the work W and black hole 
amplitude of motion and mass function A4. Keeping only 
the radius-dependent terms we find 



W (X — 



(32) 



so that the ratio W / (uorhth) is roughly homogeneous in 
1/it^ whenever <^ 1. This limit would allow to retrieve 
a scaled version of any calculation following a redefinition 
of the black hole mass and/or amplitude of motion by a 
suitable rescaling of the lengths. However the limit u ^ 1 
implies that the orbit is well inside the black hole radius of 
influence. Only a very small fraction of orbits will be found 



there. Nevertheless, (321 suggests that two configurations 



with Uofhith kept constant would yield a similar net work 
W on some orbits and so possibly the same or comparable 
imprint on the stellar energy d.f. 

We investigated this with two configurations, C2s and 
C3s, tailored to shadow runs C2 and C3 (cf. Table [l]). The 
two 'shadow' runs both had a black hole mass equal to half 
that of C2 and C3, however twice the amplitude of motion 
Uo- Table [3] lists the root mean square deviations of the en- 
ergy d.f. compared to the analytic curve (28 1. The impact 



on the energy distribution function is clearly stronger for 
the large-amplitude black hole runs. The energy deviations 
of 24% are nearly as large for an ihbh = Uo = 0.15 config- 
uration as those obtained with the same configuration but 
with rhfyh = 0.30, for which we found 29% deviations on the 
mean. In practice, we find that the energy fluctuations scale 



almost linearly with the black hole's orbital radius Uo, a con- 
clusion which we reached by comparing the results for cases 
CI to C5. What is more, the amplitude of the response of 
the stars appears robust to details of the velocity held, as 
seen when comparing conflgurations with warm or hot stel- 
lar velocity fields (Wl, HI and W2). A configuration with 
a black hole set on a circular orbit did not produce signif- 
icantly different results compared to similar configurations 
with the same mass and amplitude Uo (case Wlc, Table |3|. 
Taken together, these results highlight the importance of the 
black hole effective orbital cross section cx in transferring 
binding energy to the stars. 



5 Discussion 

Oscillations of a massive black hole about the centre of a host 
galaxy leave a signature on the kinematics of surrounding 
stars. We showed from an analytic harmonic potential model 
that stars on circular orbits lose or gain energy according to 
the relative phase between their and the black hole's orbit: 
for an evenly sampled distribution function, half of the stars 
gain energy, the other half loses energy. The black hole plays 
the role of a catalyst by allowing energy exchanges between 
the stars. 

We explored a range of orbits with black hole and stars 
in a logarithmic potential with an Bulirsch-Stoer numerical 
integration scheme (Press et al. 1992) using on the order of 
10^ orbits. The feedback of stars on the black hole orbit was 
neglected, an approach motivated by the large black hole 
to stellar mass ratio. With stars set on circular orbits, we 
found a strong response to black hole motion, at distances 
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Figure 14. Maps of the Toomre number Qj for the reference calculation C3 at A) t = 15; B) t = 16 and C) t = 17 time units of 
evolution. The shade indicates stability to fragmentation modes. The white regions have Qj < 1 and would be susceptible to fragment. 



ranging up to 3 or 4 times the black hole's radius of influ- 
ence (see Fig.[5]and[9|. Such a strong response was also seen 
in perturbed circular orbits, with velocity perturbations in 
the range ±10%. We quantified the impact of black hole mo- 
tion on the stars' energy distribution function. We measured 
root-mean square deviations growing linearly with the am- 
plitude of motion, Uo — Ro/Rc, where Rc is the core-radius 
of the logarithmic potentials. We obtained a significant re- 
sponse (signal-to-noise ratio > 4) even for modest ampli- 
tudes of ~ rbh/3, where Vbh is the black hole's radius of 
influence (cf. Eq.[5]and Table [3|. Worried that these results 
concerned circular or near-circular orbits only, we reset the 
velocity field ('warm' and 'hot' configurations, see §4) but 
kept the original energy distribution function unchanged. 
This resulted in washed out features in maps of the sur- 



face density and velocity dispersion (see Fig. 12 1 but left a 



signature on the energy d.f. of the same magnitude, with 
once again the black hole amplitude of motion the main 
agent responsible for re-shaping the energy d.f. (Table |3|. 
An analysis of the angular momentum d.f. showed that non- 
circular orbits visit the centre more frequently and couple 
more strongly with the black hole at some point on their 
orbit. Thus, although the energy exchange mechanism iden- 
tified for circular orbits plays a secondary role in systems 
with warm and hot velocity fields, the black hole motion still 
induced strong anisotropy. This was most clearly seen when 
profiling the net line-of-sight velocity v-^d resulting from two 



different projection angles (see Fig. 12 right-hand panels 



Any level of anisotropy is attributable to the motion of the 
black hole, since a static hole would have left the initial 
isotropic d.f. unchanged. 



5.1 Comparison with MW data 

The largest values of Vid were obtained from a viewing angle 
parallel to the motion of the black hole. Contrasting these 
values to the root mean square velocity dispersion, we find 
a ratio of < \vid\ > I \J < 'v\d > ~ 25% at maximum value, 
inside the hole's radius of influence (see Fig 



11 



and 



121 



Applying this to Milky Way data, where the mean velocity 
dispersion rises to ~ 180 km/s inside 1 pc of Sgr A* (Genzel 
et al.l996) we obtain streaming velocities in the range ~ 40 



km/s, a rough match to the values reported recently by Reid 
et al. (2006). The surface density proflle shows a break at 
radius ri,r ~ 0.2 pc (Schodel et al. 2007). Inside r^r, the 
volume density is fltted with a power-law index 7 ~ 1.2 
which falls outside the range 3/2 to 7/4 of the Bahcall-Wolf 
solution. Black hole motion of an amplitude Ro ~ rtr might 
cause such a break. The ratio rtr/rth ~ 0.2 compares well 
with the value ~ 0.3 adopted for our reference calculation 
(Table [1]). 



5.2 Circular black hole orbit 

We worried that a black hole set on a radial orbit might 
trigger only a subset of resonant modes from the stars, in 
contrast to the more probable situation where the hole's 
orbit has a finite angular momentum. Recall the analysis 
of §2.2, where the response was stronger for aligned orbital 
angular momenta. To test this idea, we re-ran the 'warm' Wl 
calculation with the black hole now set on a circular path 
at a radius Uo = Ro/Rc ~ 0.097, or ~ 50% its radius of 
influence (cf. Wlc, Table [l] and |3| . As the black hole orbits 
the centre, an m = density mode develops which shows 
up as a trailing arm on Fig. |12[ bottom row, middle panel. 
The black hole orbit is anti-clockwise. The arm stretches 
radially from 1 to « 2 times r^h. Its integrated mass ~ 
40% the black hole's mass, and so if the gravity of the arm 
were taken into account, the torque that this would produce 
would modify the black hole's orbit significantly, an effect 
which was neglected here. The averaged line-of-sight velocity 
of the stars, on the other hand, showed spatial variations of 
the same amplitude as in the other cases with a strictly 
radial black hole orbit of a similar amplitude. This result 
comforts the thought that black hole motion may yet give 
rise to an observable kinematic signature (especially in the 
profile of vid), regardless of the precise parameters of its 
orbit. 



5.3 Jeans instability 

Our approach suffers from a severe limitation, in that it does 
not integrate the full response of the stars to their own den- 
sity enhancements. These could become bound structures 
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which would alter the dynamics globally. To inspect whether 
this could have an influence over the evolution of the velocity 
field, we computed the Toomre number 



an 



a 

GEdl 



at each mesh points of the simulations space. Here the mesh 
size dl ~ 0.02, and the velocity dispersion a given by (311 



is measured with respect to the initial circular flow. The 
surface density E is calculated as before with an CIC algo- 
rithm. Stars are stable against self-gravitating local modes of 
fragmentation when Qj > 1.7 (see e.g. Binney & Tremaine 
1987). We applied a modifled criterion for stability, because 
the disc is presumed stable against such modes initially, that 
is, when the black hole is fixed at the centre of coordinates 
and the system is axially symmetric. We subtracted from the 
local mean square velocity dispersion the value computed for 
the symmetric configuration. In this way we measure only 
the increase in dispersion due to black hole motion, and 
set a conservative threshold for stability such that Qj > 1. 
When that condition is satisfied, the black hole contributes 
more than 58% of the square velocity dispersion required 
to prevent local self-gravitating fragmentation modes from 
growing through its orbital motion alone. Since the black 
hole already contributes more than 50% to the gravity every- 
where inside its radius of influence, it also provides the extra 
dispersion required to kill off all self-gravitating modes. 

Fig. [14] maps out Qj at three different times for the ref- 
erence calculation C3; the shaded area is for Qj > 1 with an 
upper cutoff at 2, so white means instability on that figure. 
These images should be compared to their counterpart on 
Fig. [9] It is notable that most of the over-dense thin struc- 
tures on that figure appear unstable against fragmentation 
on Fig. [m] (since they disappear in a sea of white). The outer 
dark ring at r ~ 1.3 on the figure matches the position of the 
2:5 resonance shown on Fig. [s] (using Rc — 0.46 to revert to 
physical dimensions). Thus it is very likely that structures 
that would cross this area would be heated up and disrupted 
as a result of black hole motion. This may have consequences 
for the streams of stars observed at the centre of the MW 
(Genzel et al. 2003). The dimensions of this ring, of some 
Sriih, would correspond to a radius of (roughly) 3 pc for the 
MW. This should be an element to incorporate into future 
modelling of the MW centre since actual resolution power 
already resolves sub-parsec scales. 



6 Conclusions and future work 

The response of stellar orbits is in direct proportion to the 
amplitude of motion of a massive black hole. The imprint 
of black hole motion on the stellar kinematics is in direct 
relation to the stars' angular momentum distribution func- 
tion. Stars on low-angular momentum orbit likely will col- 
lide with the hole, while those of large momenta experience 
strong beat-frequency resonances (when the hole's orbit is 
either radial or circular). The combined effect left the ve- 
locity field significantly anisotropic with a ratio of averaged 
one-dimensional velocity to rms dispersion reaching ~ 18%. 
Because analysis suggests that the black hole energy is pre- 
served while that of the stars varies in time, we say the the 



black hole is a catalyst for evolution of the stellar energy 
d.f.. 

The two-dimensional modelling done in this paper is 
a first attempt at isolating the generic features of a time- 
evolving dense nuclei with black hole motion. The quanti- 
tative outcome of the calculations would be improved in 
a study of a family of anisotropic distribution functions, 
such as e.g. the Osipkov-Merritt d.f. f{E — /rl) (Bin- 
ney & Tremaine 1987, §4.4.4), using self-consistent three- 
dimensional integrators. We have shown that when stars are 
on circular or near-circular orbits, the resonances induced by 
the black hole likely will lead to self-gravitating substruc- 
tures inside a volume of a few times rhh in diameter. Such a 
study carried out with an N-body technique is possible pro- 
vided that coUisional physics around the black hole is well 
resolved (Preto et al. 2004). Merritt (2005) and Merritt et al. 
(2007) have shown that repeated collisions with stars inside 
~ r(,h/2 lead to random walk and an effective 'Brownian' ve- 
locity transferred to the black hole. If the random walk was 
of amplitude ri,h/2, this would equally imprint the kinetic 
motion of stars outside ~ r^jh, as we have seen, through the 
catalytic process that we have outline. 

We have neglected the orbital evolution of the black 
hole. In reality the stars inside the black hole's radius of 
infiuence rbh would take away energy and lead to it sink- 
ing to the centre through dynamical friction. This does not 
invalidate the impact of black hole motion on the stellar 
kinematics because 1) this signature is manifest well out- 
side ~ 2rhh and 2) dynamical friction will be effective on a 
time-scale of ~ few orbital revolutions. Recall that the effect 
discussed here is effective on a single black hole period. 

Black hole orbital evolution would bring a higher degree 
of realism and a more fiducial comparison to observational 
data. We examined the case of a black hole on a circular 
orbit which gave rise to an m = density wave, spanning a 
mass of ~ 40% the mass of the black hole. The gravitational 
torque of the wave would rapidly brake the black hole, which 
would sink toward the centre and lock many stars along with 
it. If the wave were unstable to collapsing on itself and form a 
bound object, a double nucleus would form. The separation 
between the two nuclei would be ~ 2rbh or larger, as deduced 
from the density map on Fig. |12[ panels at the bottom. On 
the contrary, if the tidal fleld of the BH were too strong, the 
wave would merge with the black hole. This would leave the 
black hole near the galactic centre surrounded by a pool of 
stars on eccentric orbits. Tremaine (1995) has argued that 
the double nucleus of M31 (Lauer et al. 1993) may be such 
a case of an off-centre supermassive BH surrounded by a 
stretched Keplerian disc of size ~ rbh/'i (0.5" separation at 
800 pc, with Mbh ~ 8 X lO^M©). Our calculations did not 
include self-gravity, and hence the fate of the m = Q den- 
sity wave seen on Fig. [12] remains undetermined. (See also 
Peiris & Tremaine 2003, Salow & Statler 2001, and Bender 
et al. 2005 for further data on M31.) Other double or mul- 
tiple nuclei detected in external galaxies (e.g., NGC4486B 
and NGC4382, Lauer et al. 1996, 2005; the Virgo Cluster 
dwarf VCC 128, Debattista et al. 2006) are prime examples 
of the strong orbital coupling of stars with a supermassive 
black hole and its influence on the small-scale morphology 
of a galaxy. New data may reveal cases where double-nuclei 
galaxies result from the orbital coupling we have discussed 
here. For the Milky Way, current and future high-precision 
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astrometric missions, such as RAVE (accuracy of ~ 1 km/s 
out to 8 kpc at a magnitude limit of ~ 13 [/-Band] ) or GAIA 
(launch date 2011) should pick up any systematic trends in 
stellar kinematics and set firm constraints on any black hole 
motion inferred from stellar streams. 
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